function fv = rate(pro,z0,c)
% input: parameters to be calibrated including pro, z0, c
% output: equilibrium loan rate 
global beta mu_z sig_z lambda I alp m_bar eta z_min z_max rho 

%% comment out when running
% pro=pro_ini; % firm's productivity
% z0=z0_ini;

%%
R = 0.07;

z = linspace(z_min,z_max,I)';
dz = z(2) - z(1);
dz2 = dz^2;
i_entry = round((z0 - z_min)/dz);
z0 = z(i_entry);
psi = zeros(I,1);
psi(i_entry) = 1/dz;

iter = 0;
err = 1;
crit = 1e-6;
maxit = 1e3;
the = 0.01; % relaxation parameter on interest rate loop

while (err>crit) && (iter<=maxit)
    %% piecewise profit functions
    %zl = q_bar - beta*R + 1;

    qstar = z + beta*R - 1;
    piz = exp(qstar)/beta;
    
    %% HJB
    mu = mu_z*ones(I,1);
    sig = sig_z*ones(I,1);
    sig2 = sig.^2;
    
    X = -min(mu,0)/dz + sig2/(2*dz2);
    Y = -max(mu,0)/dz + min(mu,0)/dz - sig2/dz2;
    Z = max(mu,0)/dz + sig2/(2*dz2);
    A = spdiags(Y,0,I,I) + spdiags(X(2:I),-1,I,I) + spdiags([0;Z(1:I-1)],1,I,I);
    
    A(1,1) = Y(1) + X(1); % Reflecting barrier v'(z_min) = 0
    A(I,I) = Y(I) + Z(I); % Reflecting barrier v'(z_max) = 0
    
    B = (rho+lambda)*speye(I) - A;
    v = B\piz; % value functions derived from HJB
    
    %% KFE
    % entry
    m = m_bar*exp(eta*(sum(v.*psi.*dz)-c));
    %m = m_bar;

    Aa = A';
    Aa(1,1) = 0; % Reflecting barrier f(z_min) = 0
    Aa(I,I) = 0; % Reflecting barrier f(z_max) = 0
    Bb = lambda*speye(I) - Aa;
    f = Bb\(m.*psi); % stationary distribution of z derived from KFE
    %f = f/sum(f);
    K = max(sum(exp(qstar).*f.*dz),0.1 );
    
    RR = pro*alp*K^(alp - 1); % lending market clearing condition
    R = (1 - the)*R + the*RR;
    err = abs(RR-R);
    iter = iter + 1;
    fprintf('Iteration number = %4d; R = %.6f; Error = %3.0e \n',iter,R,err);
end
fv = R;

